#pragma once
extern "C" {
#include <f2c.h>
#include <clapack.h>
}
#include "LatoolException.h"
#include "IDoubleMatrix.h"

namespace LatoolNet {
	using namespace System;

	ref class DoubleGeneralBandMatrix : public IDoubleMatrix {
	private:
		long * ipiv;
		double * ab;
		int m_rownum;
		int m_colnum;
		int m_ku;
		int m_kl;
		int m_ldab;
		bool m_isFactorized;

		DoubleGeneralBandMatrix(DoubleGeneralBandMatrix ^ orig) {
			m_rownum = orig->RowNum;
			m_colnum = orig->ColNum;
			m_ku = orig->m_ku;
			m_kl = orig->m_kl;
			m_ldab = 2 * m_kl + m_ku + 1;

			ab = new double[m_ldab * m_rownum];

			for (int i = 0; i < m_ldab * m_rownum; i++) {
				ab[i] = orig->ab[i];
			}
			ipiv = new long[m_rownum];
			for (int i = 0; i < m_rownum; i++) {
				ipiv[i] = orig->ipiv[i];
			}
			m_isFactorized = orig->IsFactorized;

		};
	public:
		DoubleGeneralBandMatrix(int rownum, int colnum, int subnum, int supernum) {
			m_rownum = rownum;
			m_colnum = colnum;
			m_ku = supernum;
			m_kl = subnum;
			m_ldab = 2 * m_kl + m_ku + 1;

			ab = new double[m_ldab * m_rownum];

			for (int i = 0; i < m_ldab * m_rownum; i++) {
				ab[i] = 0.0;
			}
			ipiv = new long[m_rownum];
			m_isFactorized = false;
		};

		~DoubleGeneralBandMatrix(){
			delete[] ab;
			delete[] ipiv;
		};

		virtual IDoubleMatrix ^ Clone() {		
			return (IDoubleMatrix^) gcnew DoubleGeneralBandMatrix(this);
		};

		virtual property int RowNum {
			int get() {
				return m_rownum;
			}
		};
		virtual property int ColNum {
			int get() {
				return m_colnum;
			}
		};
		virtual property double default[int, int] {
			double get(int i, int j) {
				if ((i == j) || ((j > i) && (j - i <= m_ku)) || ((i > j) && (i - j <= m_kl))) {
					return ab[j * m_ldab + (m_kl + m_ku + i - j)];
				} else {
					return 0.0;
				}
			}
			void set(int i, int j, double value) {

				if ((i == j) || ((j > i) && (j - i <= m_ku)) || ((i > j) && (i - j <= m_kl))) {
					ab[j * m_ldab + (m_kl + m_ku + i - j)] = value;
					m_isFactorized = false;
				} else {
					if (value == 0.0) {
						return;
					} else {
						throw gcnew LatoolException("Value " + value.ToString() +  "doesn't fit to (" + i.ToString() + ", " + j.ToString() + ")");	
					}
				}
			}
		};
		virtual property MatrixType Type {
			MatrixType get() { return MatrixType::DoubleGeneralBand; }
		};
		virtual property bool IsFactorized {
			bool get() { return m_isFactorized; }
		}

		virtual void Invert() {
			throw gcnew LatoolException("Inversion of banded matrix is not supported by LAPACK.");
		};

		virtual void Factorize() {

			long m = m_rownum;
			long n = m_colnum;
			long kl = m_kl;
			long ku = m_ku;
			long ldab = m_ldab;

			long info;

			dgbtrf_(&m, &n, &kl, &ku, ab, &ldab, ipiv, &info);

			m_isFactorized = true;
		};

		virtual void Solve(IMatrix ^ b) {

			DoubleGeneralMatrix ^ gb = (DoubleGeneralMatrix^) b;

			char trans = 'N';
			long n = m_rownum;
			long kl = m_kl;
			long ku = m_ku;
			long nrhs = gb->ColNum;
			long ldab = m_ldab;
			long ldb = n;
			long info = 0;

			dgbtrs_(&trans, &n, &kl, &ku, &nrhs, ab, &ldab, ipiv, gb->dat, &ldb, &info);
		};

		virtual void SolveWithSimpleDriver(IMatrix ^ b) {

			DoubleGeneralMatrix ^ gb = (DoubleGeneralMatrix^) b;

			long n = m_rownum;
			long nrhs = gb->ColNum;
			long ldb = n;
			long info = 0;
			long kl = m_kl;
			long ku = m_ku;
			long ldab = m_ldab;

			dgbsv_(&n, &kl, &ku, &nrhs, ab, &ldab, ipiv, gb->dat, &ldb, &info);

			m_isFactorized = true;

		};
	};
}